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ABSTRACT 

Understanding the molecular content of galaxies is a critical problem in star formation and galactic 
evolution. Here we present a new method, based on a Stromgren-type analysis, to calculate the amount 
of H I that surrounds a molecular cloud irradiated by an isotropic radiation field. We consider both 
planar and spherical clouds, and H2 formation either in the gas phase or catalyzed by dust grains. 
Under the assumption that the transition from atomic to molecular gas is sharp, our method gives 
the solution without any reference to the photodissociation cross section. We test our results for the 
planar case against those of a FDR code, and find typical accuracies of about 10%. Our results are 
also consistent with the scaling relations found in Paper I of this series, but they apply to a wider 
range of physical conditions. We present simple, accurate analytic fits to our results that are suitable 
for comparison to observations and to implementation in numerical and semi-analytic models. 
Subject headings: ISM: clouds — ISM: molecules — molecular processes — radiative transfer — stars: 
formation 



1. INTRODUCTION 

1.1. Description of the Problem 

Stars form in molecular gas. It is therefore of cen- 
tral importance to understand how molecular gas forms 
in the interstellar medium of galaxies, which are per- 
vaded by far ultraviolet (FUV) radiation that destroys 
the molecules by photodissociation. In the diffuse inter- 
stellar medium (ISM) of galaxies like the Milky Way, this 
radiation is sufficient to maintain the gas occupying most 
of the volume in an atomic or ionized state. Sufficient 
concentrations of gas arc able to exclude the dissociating 
radiation, both by self-shielding and by dust absorption, 
and thereby become molecular; these are the molecular 
clouds. Regions in and around molecular clouds in which 
the thermal and chemical state of the gas is determined 
by FUV rad iation are termed photodisso ciation regions 
fPDRs: e.g. iHoUenbach fc TielensI [1991 . The surface 
density of atomic gas in a PDR, Ehi, sets the threshold 
for gas to become molecular and thus to be able to form 
stars. 

The study of PDRs is a vast field, but the ma- 
jority of the work that has been done is not suited 
to the problem of making general statements about 
the thicknesses o f HI layers. Most of it is highly 
nume r ical (Fedcrman "eral1[l979l:l van Dishoeck fc Blackl 
19861: iBlack fc van Dishoeckl 119871: iDraine fc Bertoldi 
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Browning et al.ll2003i) . which is useful for trying to match 



detailed properties of particular regions, but makes it 
difficult to draw general conclusions that allow extrap- 
olation beyond the particular numerical parameters se- 
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lected for a calculation. The same problem of how to 
draw general conclusions from particular results arises 
with galactic-scale simulations that model H2 forma- 
tion and diss ociation as one of many processes occurring 
in a galaxy (jRobertson fc Kravtsovl 120081 : iGnedin et al.l 
|2009() . In contrast, the analytic work preceeding this 
series of papers has been limited to the comparatively 
simple one-dimensional cases of a unidirectional beam of 
radiation hnpin ging on a slab do rame fc BertoldillT996l : 
ISternberd [19881 ) or a radially-con verging radiation field 
striking the surface of a sphere (jElmegreenl 119931 ) : the 
beam and slab case is a reasonable approximation for 
a PDR irradiated by a single star, but not for a cloud 
immersed in the interstellar radiation field produced by 
the ensemble of many stars. Since these authors assumed 
that the radiation is beamed, they were able to introduce 
a shielding function that describes the attenuation of the 
dissociating radiation with depth; this does not appear 
to be possible for isotropic incident radiation, particu- 
larly for spherical clouds. Moreover, these treatments 
did not identify the important dimensionless numbers 
that characterize the problem of a spherical cloud im- 
mersed in an isotropic radiation field. Perhaps the paper 
that is closest in spirit t o the present series of papers is 
that of iSternberd l|1988f ). who determined scaling laws 
for infrared fluorescent emission lines of H2 emitted by 
an irradiated slab. He assumed that the radiation was 
beamed, so that he could use a shielding function. 

We determined the dimensionless quantities that char- 
acterize ir radiated spherical clouds for the first time in 
iKrumholz et al. (2008., hereafter Paper I) , where we also 
gave an approximate solution for the atomic-to-molecular 
ratio of an irradiated gas cloud in terms of these num- 
bers. In this paper we extend this analysis by introducing 
a new and more accurate method for solving the equa- 
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tions of radiative transfer and photodissociation. We also 
extend our treatment to the case of primordial H2 for- 
mation, where there is no dust and H2 instead forms via 
the H~ channel. Finally, we apply our results to the de- 
temination of the molecular fraction in atomic-molecular 
complexes. 

1.2. The Stromgren Analysis 

PDRs are one example of a general class of problem 
in which one must find the chemical state of a gas cloud 
that simultaneously obeys coupled conditions of chemical 
and radiative equilibrium. The most famous example of 
this type of problem is that of determining the structure 
of an HII region, and the structure of PDRs and similar 
problems can be solved by ad apting the c lassic solution 
to the HII region problem bv iStromgrenI ()1939f ). In the 
absence of dust, the thickness of the layer of ionized gas 
created by an incident flux of ionizing photons F*^^ is 
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where the superscript * indicates photon quantities, rig is 
the density of electrons, and a is the radiative recombi- 
nation coefficient. Similarly, for a point source of ionizing 
radiation, the radius of the ionized gas (the Stromgren 
radius) is given by 



^st — 



1/3 



(2) 



Note that these results are independent of the value of 
the cross section, CTion, for absorption of the ionizing ra- 
diation, which is the process that actually produces the 
ionized gas. However, in order for this simple Stromgren 
analysis to be useful in determining i?st, it is necessary 
for the thickness of the layer in which the gas switches 
from primarily ionized to primarily atomic to be thin: 
^st ^ l/"-(HI)(Tion, which is generally satisfied for HII 
regions around OB stars ((Str omgren 1939). If this condi- 
tion is not satisfied, the Stromgren analysis still provides 
the total volume emission measure, / n^dV = Ll^^/a. 
The effects of dust absorption o n Stromgren layers an d 
spheres can be readily included (jPetrosian et al.|[l972[ ). 

The same type of analysis can be applied to the atomic 
gas in a PDR (Paper I). This is particularly useful be- 
cause, whereas photoionization in an HII region is a rela- 
tively simple continuum absorption process, photodisso- 
ciation is a l ine absorption process; s ophisticated calcu- 
lations fe.g.. iDraine fc Bertoldi|[l996l ) are needed to fol- 
low the absorption as the lines become opaque and even 
overlap. Only a fraction /diss of the line absorptions re- 
sult in dissociation of the H2 molecule; the rest result in 
the emission of fluorescence radiation, (/diss is a weak 
function of frequency and of the level populations of the 
H2 molecules; the avera ge value i s in the range 0-12— 
IDraine fc Bertoldi|[l996l — to 0.11— iBrowning et al.|[2003l . 
although it may be higher in regions within a PDR where 
a large fraction of the H2 molecules are in excited states.) 
The process that is analogous to radiative recombination 
is molecule formation, which occurs at a rate firi^TZ, 
where nn is the number density of hydrogen nuclei, 
/i = n(H I) /riH is the H I fraction, and TZ is the rate coef- 
ficient for H2 formation. In the Galaxy, H2 forms on dust 



rains and the rate coeffi cient is 7?. ~ 10~^^'^ cm'^ s~^ 
Draine fc Bertoldilll996( ). In primordial, metal- free gas, 
H2 forms via the H~ process or the three-body process, 
which can be described similarly (§3). In the absence of 
dust absorption, a point source with a photon luminos- 
ity L* in the Lyman- Werner bands between 91.2 nm and 
110.8 nm that is embedded in a medium of constant den- 
sity, riH, can create a sphere of atomic gas with a total 
number of neutral hydrogen atoms 
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In the absence of dust, the transition between the atomic 
and molecular gas is not sharp (Paper I); nonetheless, we 
can define a characteristic radius of the atomic gas by 
setting /i = 1, obtaining 



3/diss^ ' 



1/3 



(4) 



which has the same form as the Stromgren radius. Simi- 
larly, a dissociating photon flux Fj* incident on the sur- 
face of a slab of gas can maintain an HI column density 



fin-adz = 



./diss-^in 



(5) 



Setting Ji — 1 gives the characteristic size of this layer. 



Chi, slab 



/diss-Pjn 



(6) 



An important limitation of the Stromgren analysis is 
that the bandwidth of the dissociati ng radiation mus t 
be known in advance. As shown bv iSternberd ([1988), 
the effective bandwidth begins to shrink for metallicitics 
Z ^ Zq , where dust absorption is sufficiently strong that 
the individual absorption lines no longer overlap. For 
the cases of greatest interest, however, Z ^ Zq and the 
results presented here are applicable. 

1.3. Summary of Previous Papers in this Series 

In Paper I, we gave an approximate determination of 
the thickness of the HI layer in a sphere irradiated from 
the outside by isotropic radiation. For such a sphere, 
the incident flux is unknown, since the flux due to ra- 
diation striking the surface from outside the sphere can 
be compensated in part by the flux due to radiation that 
passes through the sphere. It is therefore advantageous 
to characterize the radiation fleld by the ambient value 
of the mean density of dissociating photons. If wc let 
E be the energy density of photodissociating radiation, 
then E* is the corresponding photon density; let Eg be 
the corresponding ambient value. We define the charac- 
teristic value of the thickness of the HI layer as ^Hi.siab 
with F*^^ replaced by cE^: 



„ _ /dissC£^0 



(7) 



For an opaque slab embedded in an isotropic radiation 
field, the energy density just outside the slab is half the 
ambient value and the flux is half that, or Fj*^ = jcEq; 

as a result, ^Hi.siab — j^ch in this case. The importance 



The Atomic-to-Molecular Transition in Galaxies III 



3 



of dust is measured by the characteristic dust optical 
depth, 

fdissO'dcE^ 



Tch = ?T-HCrd4h 



nuTZ 



X, 



(8) 



where x is one of the two dimensionless parameters in- 
troduced in Paper I. Photodissociation regions are often 
characterized by the ratio Gq/hh, where G'q is the ratio 
of the dissociating radiation field, £^01 to the typical value 
in th e Milky Way, which is 7.5 x 10~^ photons cm^"^ s^^ 
from iDraind ()1978l ). In terms of this ratio, x is 



X = 71 



fd,-21 \ G'o 
16.5/ ^^H 



(9) 



Physically, x/ fi is the ratio of the number of LW pho- 
tons absorbed by dust to the number absorbed by H2 
molecules. Thus, dust has a significant impact on the 
structure of the HI layer for x ~ 1 . The parameter y is 
of order unity in the Milky Way ijKrumholz et al.ll2009al . 
hereafter Paper II; see eg. [T3| below). C l umps in the 
PDRs discussed by iHoUenbach fc Tielend (|1999( ) have 
X ^ 0.3 — 10, while the interclump gas in those PDRs has 
values about 100 times larger. Observe that for the case 
in which H2 forms on dust grains, both <7d and TZ are ef- 
fectively measures of the total surface area of dust grains 
mixed into the gas; measures the area available for ab- 
sorbing photons, and TZ the area available for adsorbing 
hydrogen atoms. Thus their ratio should be independent 
of metallicity, dust-to-gas ratio, density, and most other 
quantities. For an opaque slab, we have seen that the 
incident flux F* - ^"^* — ^^-^ ^ " - —-^ 



cE^, so that %Lsiab = 4-4h and 
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keep in mind that these simple estimates are based on 
neglecting absorption of the incident FUV radiation by 
dust. For future reference, we note that this implies that 
the total HI column density in a dust-free slab (eq. [5]) is 
-^Hi,tot = x/(4o'd). Finally, we note that the ratio of £ch 
to the cloud radius, i?, is 



«ch 

R 



n^TZR 



X_ 

TR ' 



(11) 



where tr ~ nuadR, the dust optical depth associated 
with the cloud radius, is the other dimensionless param- 
eter in Paper I. 

As in the case of HII regions, a Stromgren-type analy- 
sis is particularly useful if the transition from atomic to 
molecular gas is sharp, £hi ^ ^/(nHC^R.^). In Paper I we 
showed that this transition is indeed relatively sharp for 
X ^ 1, so that in that case it is possible to infer the size 
of the atomic gas by setting /i 2± 1. For the case with 
dust, we make this approximation here; this is the prin- 
cipal approximation in our work. In the absence of dust, 
we can determine the total mass of H I (in the spherical 
case) or column of H I (in the slab case) without making 
this approximation. 

The advance made in Paper II was to recognize that the 
density in the atomic gas in the transition region is set 
by the req uirement that it be in two-phase equilibrium 
(jWolfire et al.„2003i) . The latter authors showed that the 



minimum density of cold HI (the CNM) in galaxies is 
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where Eq and Z' are the radiation field and metallicity 
normalized to the local Milky Way value. If the density 
is several times the minimum value (nn = (/'cnm'T'H, min 
with ipcNM 3), then 



X = 3.1 



c^d, -21 



^-16.5('/'CNm/3) 
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/0.365 



4.1 



(13) 



depends primarily on metallicity, and weakly at that; 
here ad, -21 = ad/i^O'^^ cm^), etc. The surface den- 
sity of H I in a galaxy, assuming that it is in a slab and 
neglecting dust absorption, is then 



^HI — THI.slab, 

(Td 



^'(0cnm/3) 



1 + 3.1Z 



/0.365 



4.1 



Mq pc-2 (14) 



for TZ-ie.5/(Jd,-2i = 1; here fin ~ 2.34 x lO^^'^ g is the 
mass per hydrogen for a gas of cosmic abundances. Pa- 
per II showed that dust absorption reduces this by about 
a factor of 2. On the other hand, if the slab of gas is illu- 
minated from both sides, then the total column of HI is 
increased by a factor 2. Paper II generalized these results 
to the case of finite clouds and demonstrated that the re- 
sults were in good agreement with existing observations. 
Indeed, since the theory applies to individual cloud com- 
plexes, predictions were made on how the results would 
c hange as the reso l ution of the observations improves. 

iKrumholz et al.l (|2009bl ) applied the results of Pa- 
pers I and II together with the earlier work of 
IKrumholz fc McKed (|2005[ ) to determine the star forma- 
tion rate in galaxies as a function of metallicity. For 
low surface densities (E 10/ Z' Mq pc"^), the rate is 
dominated by the transition from HI to H2; for interme- 



diate columns {10/Z' Mq pc 



S < 85 Mq pc" 



it 

is determined by the properties of giant molecular clouds; 
and for larger column densities it is de termined by the 
pressur e of the galactic ISM. Similarly, IKrumholz et al.l 
(|2009c[ ) showed how the results of Papers I and II can be 
used to explain the observed absence of damped Lyman- 
a systems with high column densities and metallicities. 
The treatment of the HI — H2 transition in both these 
papers relied on an analysis that is more broadly appli- 
cable than the one developed in Paper I. We provide that 
analysis here. 

In §2, we present the general formalism for determin- 
ing the thickness of the HI shielding layer for spher- 
ical molecular clouds. This formalism is based on a 
Stromgren-type analysis. In §3, we present the results 
for both slabs and spheres in the absence of dust; this is 
relevant to the study of primordial clouds. §4 gives the 
results including dust absorption, §5 shows how these 
results can be generalized to the case of clouds in which 
the atomic and molecular regions have differing densities, 
and the conclusions are given in §6. 

2. FORMALISM 
2.1. Basic Equations 
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Consider a spherical cloud of radius R and uniform 
density exposed to an isotropic UV radiation field. We 
describe the radiation field in terms of the specific inten- 
sity in photon units, /* = Iv/hv. Let = n2iJ2v be the 
opacity due to absorption by H2 molecules, where 712 is 
the density of H2 molecules and <T2v is the H2 absorption 
cross section at frequency v. Let = n^^ad be the opac- 
ity due to absorption by dust, where rtn is the density of 
H nuclei and ct^ is the dust absorption cross section per 
H in the photodissociation part of the spectrum, from 
91.2 nm to 110.8 nm; we ignore the weak frequency de- 
pendence of this cross section. The equation of transfer 
is then 



-{k2i, + Kd)K = -("20-2iy + nu(7d)K (15) 



for the intensity in the e direction. The opacity K2u de- 
pends on position, but we have suppressed the argument 
for clarity. For a cloud bathed in a uniform radiation 
field of intensity /*o, the intensity at a point inside the 
cloud is 

K = Ko cxp(-r2^ - r), (16) 

where the optical depths are proportional to the column 
density: the dust optical depth is 



ii-aCTdds 



and the H2 optical depth is 



T2u 



Jl2(J2iyds ~ N2(J2v: 



(17) 



(18) 



where the range of integration extends along the ray from 
where the ray enters the cloud to the point in question 
and N2 is the corresponding H2 column density along the 
ray. 

We next integrate over the range of frequencies that 
can photodissociate H2, i^2 ^ ^ '^ij corresponding to 
wavelengths 91.2 nm < A < 110.8 nm. The frequency- 
integrated equation of transfer is 



V/* 



-[Ki{e) + Kd]I* 



where 



K2i,I*diy. 



(19) 



(20) 



We include the argument e in kj to emphasize that this 
opacity depends on angle through its dependence on /* ; 
Kj also depends on position. The intensity-weighted H2 
opacity is then a function of N2, 



Kl 



K21' exp(-T2i.) dv 
Q exp(-r2i.) dv 



Kl{N2 



(21) 



The dependence of kj on N2 is governed by two coun- 
tervailing factors. On one hand, as N2 increases more 
and more photons at frequencies near line centers are 
absorbed. Since a2i> decreases away from line center, the 
mean the opacity per H2 molecule is reduced. On the 
other hand, K2jy is also proportional to the density of H2 
molecules n2, which increases with N2. The intensity- 
weighted optical depth that determines the attenuation 
of /* is T/ = / Kj ds, so that 



r =/*cxp[-r/+r]. 



(22) 



In terms of the first three angular moments of the ra- 
diation field, 

E; = Jdrt f; = J dn e/;, v^-J an eei;, 

. (23) 

the first two angular moments of the equation of radiative 
transfer arc 



V f: 

VP*, 



-c{k2i, + Kd)E*^ 



Integrating over frequency gives 

V ■F*^-c{KE + tid)E\ 



V P* 



-{KF + Hd)F* 
c 



(24) 
(25) 

(26) 
(27) 



where 



KE = 



dvK2yEl 



E* _ 

j d^jJdyK2j: ^^Jdn Kjl*. (28) 

We have assumed that the direction of F* is independent 
of frequency; as a result, is given by an expression like 
that for KB with E* replaced by \F* \ and with a factor jj. 
in the angular integrals. Writing out equations and 
(j27p in the case of spherical symmetry, wc obtain 

1 d 



■r'F*=c{KE + Kd)E*, (29) 
dP* 3P* - E* 1 



dr 



-{kf + Kd)F*, 



(30) 



where F* = = -F* and P* = P*^ is the rr compo- 
nent of the radiation pressure tensor. 
Now define the characteristic flux 



^ch = —7 1 

/diss 



(31) 



which is the flux required to photodissociate a layer of 
thickness i?, as discussed in §1. The solution is deter- 
mined by the two quantities defined in Paper I, 



TR = nYi(JdR, 



(32) 
(33) 



where E^ is the number density of dissociating photons 
far from the cloud. In terms of these parameters, we have 



(34) 



We define the normalized fiux and photon density by 

cE* 



E* 



(35) 



F, 



ch 



The normalized photon density far from the cloud is then 

X 



E* = ^. 
tr 



(36) 
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We now use the assumption that the formation and 
destruction of H2 are in balance in order to determine 
the opacity ke] this is the essence of the Stromgren ar- 
gument. Observe that it is difficult to calculate ke ex- 
plicitly, since it involves first a frequency integration over 
the cross section and the specific intensity, both of which 
are complicated functions of frequency due to line ab- 
sorption, and then an average over angle. Balancing the 
rates of formation and destruction of H2 gives 

dy K2„I*, (37) 
(38) 

where /i = n(HI)/nH is the atomic fraction and TZ is the 
rate coefficient for H2 formation. The principal approx- 
imation in our analysis is to assume that the transition 
from atomic to molecular gas is sharp, so that /i ~ 1 in 
the atomic gas. The balance equation then implies that 



J J ui 

^ fdissCKsE* , 



Ke 



1 



cfd\ssE* RE* 



(39) 



Our scond approximation is to relate the flux-weighted 
opacity to he by 

He = (pK^E, (40) 

where 0k is a constant. The results for the ID case below 
verify that this is a good approximation. We expect 0^ to 
exceed unity since the rays that carry the flux are closer 
to the radial direction and have therefore undergone less 
attenuation and interact with a larger cross section, as 
mentioned below equation (j2ip . 

To close the moment equations of radiative transfer, 
we introduce a variable Eddington factor. 



/Edd — 



E* 



(41) 



Expressing equations (|29p and (j30p in dimensionlcss form 
and using equation (j39p to eliminate k^;, we obtain 



dF* 2F* ^ 
dx X 



-7- (/Edd£^* 

dx 



(3/1 



Edd 



1)E* 



E* 



(42) 
trF*{AZ) 



where x = r/R. These equations must be solved subject 
to the boundary conditions that both the photon density 
and the flux vanish at aiHa : the transition point from 

atomic to molecular gas; i.e., E*{xb_2) = F*{xYi2) = 0- 
The values of the photon density and flux at the surface 
are not known, except in the case in which the cloud is 
completely opaque. 

2.2. Approximate Determination of the Eddington 
Factor /Edd 

Many astrophysical problems are characterized by a 
central source of radiation in an opaque region, where the 
radiation is nearly isotropic, surrounded by a less opaque 
region in which the radiation is beamed outwards. Shu 
(1991) writes a general form for a relation to close the 
moment equations of radiative transfer as 



E* 



3P>* _ S—E* 



Free space 




Fig. 1. — Diagram of the geometry of a spherical eloud illumi- 
nated by an external radiation field. 



He points out that the choice qe = 2 makes a smooth 
transition from the isotropic case {E* ~ 3P* , F* — > 0) 
to the beamed case (E* ^ P* ^ F* /c). Our problem 
is quite different, however: the radiation is isotropic at 
large distances from the cloud, but close to the cloud 
there is a "hole" in the radiation field since radiation 
cannot pass through the molecular core. We therefore 
need to develop a different approach to closure. 

As remarked above, we close the moment equations by 
using variable Eddington factor, /Edd- Our approxima- 
tion for /sdd at any point inside the cloud is based on 
representing the main qualititivc feature of the radiation 
field — that the intensity decreases for rays at larger an- 
gles to the outward normal and vanishes for rays that 
have intersected the molecular core — by a step function 
decrease at some angle 9f. (We adopt the convention in 
Paper I that 9 is the angle between the outward normal 
and the direction from which the photon originates - see 
Figure [TJ) We thus approximate the radiation field as 
being of constant intensity for < < 6 f (correspond- 
ing to 1 > ^ > /i/, where ^ = cos 9) and vanishing for 
6f < 9 < n (corresponding to /i/ > /.j > —1). We then 
obtain 



/Edd 



'e* 



1 



(45) 



(44) 



It remains to adopt an ansatz for fxj. First consider 
the dust-free case. For points close to the molecular core 
{x we assume that 29 f is the angle subtended 

by the molecular core, so that sin9f ~ Xfja/^;. In the 
opposite limit, in which the molecular core is assumed 
to be small (xhs ^ 1) and the point is far from the 
core, we assume that a ray at 9f intersects the midplane 
at a distance XfR from the center of the cloud, so that 
tan 0/ — Xf jx. We combine these two limits by assuming 
that tan^ 6*/ is the sum of the squares of the tangents in 
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the two limiting cases, 



tan^ 9f = 



4. 



(46) 



The corresponding cosine is — —1/(1 + tan^fl/)^/^, 
where the minus sign reflects the fact that the angle is 
in the second quadrant: 



x"^ + Xj{x'^ - xf^J 



1/2 



(47) 



The parameter xj is an eigenvalue of the problem, which 
is adjusted to give the correct net flux at the surface. 

How does this compare with Shu's (1991) closure ap- 
proximation? Using the model that the radiation field is 
uniform for fi > and zero for fi < iij implies that the 
flux is jE*{l — ^'j). Solving for the closure factor qp in 
equation (jH]), we find qp = '^Hf- Since /i/ is negative, it 
follows that the closure factor is negative, which is qual- 
itatively different than the standard case. The fact that 
< implies that /Edd in equation (|45p is always less 
than i; it is in the range i > /Edd > j- 

Next, consider the case in which dust is included. For a 
dusty slab, the dissociating photon density at an optical 
depth T from the surface is 



1 



Evaluation of the integral gives 



2"0 1^ exp(-T//i)d^i, 



(dusty slab). (48) 



E* ^ -E*E2(t) ~ -E* — 
2 / 2 " H 



(49) 



where E2 is an exponential integral and where the ap- 
proximation is accurate to about 25%. The advantage of 
the approximation is that it shows that the photon den- 
sity is as if the intensity along the normal filled a solid 
angle corresponding to Afl = 27rA/^ = 27r(l — /i^) with 
1 — /id = 1/(1 -I- r) , so that 



are attenuated, and those that strike the molecular core 
(m < M2) arc absorbed. Integration of equation (|22|) then 
gives 



E*{rn) = Ie* 



1 



dn e 



-(ti+t) 



(53) 



We then introduce our next approximation: we replace 
T/ = J Kj{r, fi)ds, with te = J KEir)ds (see eq. [28)) . 
Note that the error introduced by this approximation is 
limited, since the term in which it is made contributes at 
most half the total. Proceeding in the same fashion with 
the evaluation of the flux at the surface, we obtain 



E*{rR) 



1 + f dfi e-("^+^) 



l + 2/^%/./xe-( 



te+t) 



(54) 



where the factor 2 in front arises since the unattenuated 
intensity is ^Eq whereas the unattenuated flux is jcEq. 
The optical depth is given by 



TE 



x'dx' 



(x'^-sin^ ey/^E* 



(55) 



where sinO is the value of x of closest approach to the 
center of the cloud for a ray at an angle 6. 

The equations we have now written down constitute a 
complete set that fully determines the solution. For a 
given choice of x a-nd r^j, a solution consists of two un- 
known numbers, fif (or equivalently /Edd, since the two 
are related by equation I5ip and xn.^ , and two unknown 

functions, E* and F* . The two functions are constrained 
by the ordinary differential equations ([42]) and (j43|) , while 
the two numbers are constrained by the algebraic equa- 
tions E*{x}i^) = F*{x-ii2) — and the consistency con- 
dition equation ([54]) . We must choose a fixed value of 
4>K, since we lack an equation to determine an additional 
parameter. We discuss how to choose this value in i; 12.41 
and we defer discussion of how to obtain the solution 
numerically once we have chosen 0^ until 13. 21 and 14.21 



Aid = 



1 



(50) 



For the spherical case, we take r = (1 — x)t/j, where tr 
is the dust optical depth from the cloud center to edge. 
To join the case where dust dominates the opacity onto 
the case where molecular absorption dominates it, we use 
the simple ansatz 



/Edd = i(l -l-Ait +M?)- 



(51) 



with fj,t = fif + fid- Note that whereas /i/ is always 
negative in the dust-free case, it can be positive when 
dust is important. 

2.3. Boundary Conditions for the Spherical Case 

As in Paper I, we use the solution of the radiative 
transfer equation at the cloud surface to constrain the 
solution. At the surface {x = 1), photons at ^ > are 
unattenuated, those that traverse the cloud outside the 
molecular core, > yu > /i2, with 



1 



.1/2 



(52) 



2.4. The Semi-Infinite Slab Limit 

The normalization we have used for the spherical case, 
which is based on F*^ cx R, breaks down for a one- 
dimensional, semi-infinite slab, which corresponds to the 
limit i? — > 00. We therefore normalize with respect to 
the ambient radiation field, E^ by defining 



cEX 



El ^ — ; 
' E*' 



(56) 



i at the surface of the slab. 



thus, Fi*(0) = i and El{Q) 

The first two moments of the radiative transfer equation 
become: 



dFl 

rf(/Edd-E*) 



---El, 
X 

XEI 



F*, 



(57) 
(58) 



where dr = nuaddz, with z being zero at the surface 
and increasing inward. In the ID limit (x, xn^ — > 1), the 
angle /if that enters the Eddington factor approaches 
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zero (eq. ITT]) . The Eddington factor thus becomes 

/Edd = ^(l + /i<i + Md)- (59) 

For small r (i.e., close to the surface), /sdd ^ | as it 
should, since the radiation is approximately isotropic in 
the half-space ^ > 0. For large r, we have /x^ = t/(1 + 
r) — > 1, and /Edd 1; this too is as it should be, since 
highly extincted radiation is beamed so that P* ^ E* . 

As with the spherical case, we have now written down 
enough equations to close the system. The independent 
parameter that determines the solution is the dimension- 
less radiation intensity x- For a given x, a solution con- 
sists of two numbers, the ratio (pK. of the flux-mean opac- 
ity to the energy-mean opacity and the dust optical depth 
at which the molecular transition occurs, plus two 
functions, and F* . These quantities are constrained 
by two algebraic equations, El{T]i^) = F*(tb_2) = and 
the two ODEs ^ and ([58]). These ODEs are subject to 
the boundary conditions £^1(0) = 1/2 and F^{0) = 1/4. 
This set of constraints therefore fully specifics the sys- 
tem. We defer discussion of our solution algorithm to 
§ 13.11 and 14.11 We pause here to point out an important 
distinction between the spherical and slab cases: for the 
spherical case, the angle 0/ at which we transition from 
zero flux to finite flux represents an extra parameter to 
be determined that is not present in the slab case. As a 
result, for the slab case we can use our constraint equa- 
tions to determine 4>k , and we can then use the resulting 
value of in the spherical case. As we shall see, a single 
value of works quite well over a very broad range of 
radiation fields. 

3. THE DUST-FREE CASE: PRIMORDIAL CLOUDS 

Although the motivation for this series of papers was to 
understand the atomic-molecular transition in galaxies 
today, where II2 is formed on grains, our methodology 
applies to primordial gas clouds also, where H2 is formed 
via gas phase reactions. The effective values of TZ for 
gas-phase production of H2 are given in the Appendix. 
The results obtained here are also useful for comparison 
with the results for dusty clouds in the following section. 

3.1. Semi-Infinite Slab 

In the limit that — > 0, equations ((57|) and ([58]) for 
the ID case become 



dFl 
dNHi 

dfEddEi 

dNm 



1 

'4' 
iE* 



(60) 
(61) 



where we used the fact that Kd/x — "-HCd/x 
nH/(4A^Hi.tot) (see below cq. [T0|) and defined 



HI 



HI 



(62) 



which varies from at the surface of the slab to 1 at the 
point that the slab becomes fully molecular. Note that 
we have not assumed anything about the spatial varia- 
tion of the HI fraction, /i. As pointed out in ij2.41 the 
boundary conditions are Fi{l) = j and E^ = i; we 
also pointed out that for small t, the Eddington factor 



approaches /sdd =5- In the dust-free case, the use of 
/Edd = at large optical depths in the H2 lines is justi- 
fied by the fact that the photons that dominate the pho- 
todissociation arc in the line wings; they arc not highly 
extincted and therefore are not beamed. Integration of 
the equations then gives 



F* = -(l-iVHi), 



1 



30K 



^ 4 8 * 



mi 



HI 



(63) 
(64) 



Since E^ must vanish at the same point that F* does 
(i.e., at Nui = 1), it follows that 



(65) 



We shall see in H4.1l that this remains a good approxima- 
tion in the dusty case for all physically relevant values of 



3.2. Spherical Clouds 

The problem of the dust-free spherical cloud does not 
have an exact solution with the Stromgren method be- 
cause the net incident flux, F*^^ , depends on the structure 
of the molecular hydrogen in the transition zone; some 
of the rays that strike the cloud will penetrate all the 
way through, reducing the incident flux by an amount 
that depends on the unknown spatial distribution of the 
H2. Wc therefore consider two complementary approxi- 
mations: In the flrst method, wc assume that the gas is 
entirely atomic in the transition zone (/i = 1) and cal- 
culate the incident flux self-consistently. In the second 
method, we allow for the fact that molecular gas exists in 
the transition zone (/i < 1), but assume that the tran- 
sition zone is thin enough that the incident flux has the 
value appropriate for an opaque cloud, -F'i*^. ~ \cEq. We 
shall see that the results overlap, giving a self-consistent 
determination of the region of validity of our results. 

In general, both x and go to zero for dust-free spher- 
ical clouds, since both are proportional to the dust cross 
section ad- However, the ratio 



/d 



f = A - ^£Ji - 



issCEq 



TZR 



(66) 



remains flnite, and this becomes the independent param- 
eter describing the normalized radiation intensity in the 
dust-free case. We anticipate that the cloud will be fully 
atomic at high values of ^, whereas the atomic gas will 
be confined to a thin shell at low ^. 



With TR 



3.2.1. Method 1: fi^l 
and /i = 1, equation (|42)) reduces to 



dF* 2F* 
dx X 



1 



(67) 



which may be integrated analytically, subject to the 
boundary condition F*(xii2) — 0, to obtain 



F* 



1 



X 



(68) 
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Fig. 2. — Dust-free sphere. Fraction of volume occupied by 
molecular gas {upper panel) and ratio of molecular to atomic 

volume [xii^/{l — xn^)]^ {lower panel) as a function of 5 = xItr- 
The solid line is the numerical solution, and the thic k dashed line 
is the analytic approximation given by equation II70I I. 



Substituting this into equation (|43p gives 

— (/Edd^*) + ~ = "^"^ - ('^hJx)3] 

dx X E* 

(69) 

We can now obtain the relationship between xn^ and 
^ via the foUowing algorithm: 

1. Select a value for XHa- Our algorithm will deter- 
mine the corresponding ^. 

2. Pick a trial value of the parameter Xf. (The value 
of is fixed to 4/3, as discussed in § 13.11 ) The 
choice oi Xf determines the Eddington factor /sdd 
via equation ((45|) . 



3. Given /Eddi determine the function E* by inte- 
grating the ordinary differential equation (|69p from 
a; = .th2 to .t = 1, using the boundary condition 
that E*[xu^) = 0. 

4. Check whether the solution for E* obeys the inte- 
gral constraint equation (|54p to within some spec- 
ified tolerance. If not, repeat steps 2-4, using 
Newton- Raphson iteration to find a value oi Xf 
that minimizes the difference between the two sides 
of equation ([5^ . 

5. If the solution does satisfy the constraint equation 
dSlD, equation ^ tells us that E*{1). 

Applying this algorithm yields the curve for xb_^ versus ^ 
shown in Figure [21 Note that in this approximation, 
goes to zero exactly for ^ = 2. In an exact calculation 
with /i < 1, the gas would not be fully atomic at this 



point, but this result shows that clouds with ^ > 2 are 
substantially atomic. The formula 



1 



1 + 0.0712^2.8 



(70) 



matches the numerical solution to within 11% in the 
range > 0.01, and may be used in place of a nu- 
merical evaluation for most practical calculations. 



3.2.2. Method 2: F-* 



If we retain the dependence on /i but assume that the 
incident flux is known, an analytic solution is possible. 
In this case, equation (|^^ reduces to 



dF* 2F* 
dx X 
The boundary condition F^^^ = 



F*(l) 



(71) 



jcE^ becomes 



(72) 



with the aid of equation (j36p . Integration of equation 
dZH) gives 



x'^fidx = 



A/h 



= 3 ( ^ - x^F' 

loud V 4 



(73) 



where A/hi is the number of HI atoms outside r — xR 
and A/h, cloud = 47ri?'^nH/3 is the total number of H nuclei 
in the cloud. To relate this result to that from the first 
method, we define an effective value of xn^ such that 

A/h 



3 , -'Vlli,t0t /r,.\ 

xii2,es = ^- J? , (74) 

A/h, cloud 

where Ah I, tot J the total number of HI atoms, is evaluated 
at a;^F*/(^/4) ^ 0. Equation ^ then imphes 



^H,.cff 



(75) 



This has the same form as equation (fZO)) . which is based 
on the first method, for ^ ;S 1- Since the first method 
properly accounts for the flux that penetrates through 
the cloud, we conclude that the total number of HI atoms 
in a spherical, dust-free cloud is 



A/h I, tot — ^^A/h, cloud 



(76) 



for ^ 1. Since Method 1 showed that the cloud be- 
comes fully atomic for ^ ^ 2 and both methods agree for 
^ ^ 1, we conclude that the results of this section are 
valid everywhere except where the cloud is close to being 
fully atomic. 

4. DUSTY CLOUDS (Fi = 1) 

4.1. Semi-Infinite Slab 

We now consider the case of clouds with dust, (Jd > 0, 
beginning with the semi-infinite slab. As discussed in 
§1, the presence of dust sharpens the transition between 
atomic and molecular gas, enabling us to make the ap- 
proximation that /i = 1 in the atomic gas. Our goal is to 
determine the relationship between the imposed dimen- 
sionless radiation field x and the dust optical depth THa 
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Fig. 3. — Upper panel: dust optical depth of the atomic-to- 
molecular transition th^ as a function of radiation field x- The 
solid line is the numerical solution and the thi ck dashed line is the 
analytic approximation given by equation I I77I I. Symbols show rjjj 
versus x for the Meudon runs, with circles showing the results if 
we adopt the integral definition of th^ and Xi and triangles for the 
point definition. See i) l4~l.ll for details. Lower panel: log{0K — 4/3) 
as a function of x- Paper II suggested that galaxies typically have 
X ~ 2 -- 3 for metallicities Z ~ (0.1 — ^)Zq. 

at which the transition from atomic to molecular occurs. 
In the process we must determine the unknown param- 
eter (f>f^. We do so using an algorithm slightly modified 
from that in 



1. Select a value for x- Our algorithm will determine 
the corresponding ■ 

2. Pick a trial value of the parameter . 

3. Given (j)^^ integrate the ordinary differential equa- 
tions ((57)) and ([55]) using the boundary conditions 
that El{0) = 1/2 and F^{0) = 1/4. Integrate 
in the direction of increasing r, stopping once we 
reach a value r = tq for which i<\*(ro) = or 
-El* (to) =0. 

4. If F^{tq) = 0, check whether EKtq) = as weh 
to within some specified tolerance; if EKtq) = 0, 
check whether F^{to) = to within the specified 
tolerance. If not, repeat steps 2-4 using Newton- 
Raphson iteration on to minimize the value of 
EtiTo) orFi*(To). 



5. If Eliro) = Fi*(to) 
tolerance, then th, = 



to within the specified 



To- 



Applying this algorithm produces the curve for as a 
function of x shown in Figure [3) The formula 



TH2 



In 



1 



X 



(77) 



4(l + e)_ 

with e = 0.122 reproduces the numerical solution 
to within 2% for logx < 2.5. 



It is interestin g to compare our result with that of 
iSternberS ()1988[ ) , who wrote the chemical balance equa- 
tion (1381) as 



/ln|7^ = i?/.hicld(iV2)n2e-('^^+''^^'"^ 



(78) 



where Dn2{0) = /dissCK£;(0)i?*(0) is the dissociation 
rate at the clo ud surface (this is a slight chan ge from 
Sternber"^ll988l . who defined D as the rate for the lDrain^ 
19781 radiation field). Integration of this equation, with- 
out the assumption that /i ~ 1, gives 



a,/^ n,dz^ In , ^^^^ 



DG 



(79) 



where ni is the atomic hydrogen number density and 



G = ad 



/shicid(Af2)e-2^^-^dAr2 



(80) 



Provided the transition from atomic to molecular gas is 
sharp, which is a good approximation in the dusty case 
when X 1 (Paper I), then the LHS of equation (|7^ 
is just THa- Sternberg's result then has the same form 
as ours, although it requires knowledge of the shielding 
function /shield hi order to evaluate it accurately. Fur- 
thermore, Sternberg considered the case of beamed ra- 
diation, whereas in our problem the ambient radiation 
is isotropic far from the cloud. For beamed radiation, 
F* = cE* , so that equation ([77]) becomes — In(l-l-x) 
(neglecting e); Sternberg's factor DG/uhTZ is thus com- 
parable to our factor Xi and it is possible to show this 
directly (Sternberg, private communication). In Paper 
I we compared the results from our Stromgren analy- 
sis for beamed radiation with those from a detailed nu- 
merical calculation of H2 photodissociation and found 
good agreement. Our results are therefore consistent 
with those of Sternberg. 

It is also interesting to examine the value of deter- 
mined by our procedure. For the dust-free case 0^ = 4/3 
exactly, and the lower panel of Figure [3] shows that this 
remains a very reasonable approximation even when x 
is significant, meaning that dust opacity plays an im- 
portant role in shielding molecules. The approximation 
(/)^ = 4/3 is good to better than 10% for x < 12.6, and 
good to within 50% for x < 50. Since typical molecular 
clouds have x — 1 (Paper II), this means that we can 
adopt = 4/3 in general without significant error for 
the parameters that describe real clouds. 

4.1.1. Comparison to Numerical Results 

The dusty slab case is a useful one to consider be- 
cause there are a number of PDR codes that solve for 
the transition for HI to H2 in an infinite slab exposed 
to an isotropic external radiation field using brute-force 
frequency-dependent radiative transfer through the LW 
bands, coupled to explicit calculation of the H2 level pop- 
ulations. These codes also include a number of physical 
processes that our treatment omits, such as HI produc- 
tion by cosmic rays, temperature dependence of the H2 
formation rate coefficient, and variations in density and 
temperature through the PDR. By checking our simple 
analytic model against these much more sophisticated 
(and computationally costly) calculations, we can esti- 
mate the likely level of error in our method. 
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For this comparison we examine a series of calcula- 
tions performed by M. G. Wolfire (2009, p rivate commu- 
nicat i on) using the Meudo n PDR code (|Le Petit et alj 
120061 : iKaufman et al.ir2006n . In these models the radia- 
tion field outside the slab is set to Gq = 0.3, 0.5, 1.5, 
3, 10, 30, or 100, where Gq is measured in units of the 
iDraind (|1978[ ) field. The dust cross section in the Meudon 
runs is (Td = 1.1 X 10~^^ cm^. The H2 formation rate co- 
efficient varies with temperature, but a rough mean value 
is TZ = 3 X 10~^^ cm^ s~^. Most runs use a cosmic ray 
ionization rate ( = 10^^^ s^^ (the canonical Milky Way 
value), but one also uses C = 10^^^ s~^, a factor of 10 
higher. The gas density in the Meudon code can be set 
to a constant value, but to provide a more stringent test 
for how our method fares when there is a realistic level of 
density variation, we instead compare to runs where the 
density is determined in one of two self-consistent ways. 
In some runs the density is determined by requiring that 
the thermal pressure in every cell have a constant value 
Pth/ks = 10^ or 10^ K cm~'^; the density and tempera- 
ture are then obtained iteratively by requiring that the 
gas be in mechanical and thermal equilibrium at this 
pressure. In the second group of runs, the mean density 
as a function of depth in the molecular cloud is deter- 
mined by requiring that the gas be in hydrostatic equi- 
librium with the self-gravity of a molecular cloud of a 
specified mass Aid- In these models the median density 
also differs from the mean density due to the effects of 
turbulence-induced clumping. For these models the den- 
sity at a given depth is fully determined by the cloud 
mass Mci, which varies 10^ Mq to 3 x 10*^ Mq. Full de- 
tails of the models are given in Wolfire, HoUenbach, & 
McKee (2009, in preparation). 

In order to compare the numerical results to our ana- 
lytic predictions, wc must choose definitions for ambigu- 
ous quantities. First, our model is for constant-density 
gas, but the density is not constant in the numerical runs. 
Thus we must choose some sort of average density to 
plug into our analytic formula for x- This choice is non- 
trivial because the density rises sharply as gas converts 
from H I to H2 in the numerical runs, increasing by an or- 
der of magnitude between the predominantly atomic and 
predominantly molecular region. Clearly our predictions 
will depend on whether we take the gas density to be the 
atomic density, the molecular density, or some interme- 
diate value. Second, the atomic-to-molecular transition 
is fairly sharp, but the atomic fraction never falls to zero 
exactly in the numerical calculations. Thus we must de- 
cide how to define the column density or optical depth 
of the atomic layer. 

We consider two schemes for how to define these two 
quantities. First, following our approach in Paper I, we 
can use an integral definition of the optical depth of the 
shielding layer: 

TH2 ^ cTd J ni dz, (81) 

where ni is the number density of atomic hydrogen. In 
the limit of an infinitely sharp transition from atomic to 
molecular, this clearly gives the correct optical depth for 
the HI layer, and this definition has the advantage that 
we need not choose an arbitrary molecular fraction at 
which to declare that the transition has occurred. Since 



the HI fraction falls off faster than exponentially, the 
integral is guaranteed to converge (neglecting the tiny 
amount of HI produced by cosmic rays). For this defi- 
nition of , since we are sampling the entire gas col- 
umn where the atomic to molecular transition occurs, 
we choose to define the density to be the mass-averaged 
density of all computational cells with < 0.9, where 
= 2n2 /nn is the fraction of the gas in the molecular 
phase. Thus our definition takes the density to be the av- 
erage density considering both the predominantly atomic 
gas and the gas in the transition zone where the molec- 
ular fraction is below 90%. However, we also considered 
alternate methods to compute the mean density, ranging 
from considering only cells with /h^ < 0.5 (i.e. exclud- 
ing more of the transition zone) up to allowing all cells 
with < 0.95 (i.e. going further into the molecule- 
dominated region). These alternate definitions do not 
change the predicted value of by more than a few 
percent. 

The second option we consider is to define as the 
dust optical depth from the edge of the slab up to the 
point where the /h2 = 0.5. We refer to this as the point 
definition, since we are choosing to measure the column 
up to some specific point. The choice of 0.5 is rather ar- 
bitrary, but picking a particular point allows us to study 
how big an error we might make by assuming that the 
atomic-to-molecular transition is sharp, something that 
an integral definition such as equation (jSip might ob- 
scure. For this point definition of th^ , we need to define 
the density differently than we did using the integral def- 
inition - clearly the density of gas in which is in the 
range 0.5 — 0.9, which we include in our averaging for the 
integral definition, cannot affect the propagation of radi- 
ation in the region where fn^ < 0.5. Instead, since the 
point where /hj ~ 0.5 is primarily determined by how 
much radiation is absorbed in the predominantly atomic 
region, we estimate the density in this case by computing 
the mass-averaged density of gas in which fn^ < 0.05 (i.e. 
in the mainly atomic region). If we instead consider gas 
with /h2 < 0.025 or gas with f^^ < 0.1, our predicted 
optical depths th^ change by ^ 50%. 

Now that we have defined a mean density tt-h , we must 
compute the other quantities that determine x- the ra- 
diation field Eq and the H2 formation rate coefficient TZ. 
For the radiation field, the photon number density Eq 
corresponding to a radiation intensity Gq at the cloud 
surface is E^ = 2G'q{7.5 x lO""^ cm^^)^ where 7.5 x 10"^ 
cm~fis the LW photon density corresponding to the 
IDraind (|1978[ ) field, and the factor of 2 arises because 

the free-space value of the radiation field, Eq, is twice 
the value at the surface of an infinite slab. For the H2 
formation rate coefficient, we use 7^ = 3 x 10~^^ cm'^ 
s~^, the rough mean value given the temperature varia- 
tion within the PDR. Finally, we have ad = 1.1 x 10~^^ 
cm^ and /diss = 0.11. Given these values, we compute 
X for each Meudon run, and then we compute the dust 
optical depth of the H I - H2 transition predicted by our 
analytic model using equation ([77)1 . 

The circles in the upper panel of Figure [3] show TH2 
versus x for integral definition for and density, 
while the triangles show the results if we instead adopt 
the point definition. In the range x 1, the choice of 
definition matters little and the errors are small. In this 
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region the largest absolute difference between our ana- 
lytic model and any of the numerical results is 23% or 
28% and the median absolute differences are 4% and 9% 
for the integral and point definitions, respectively. Over- 
all our ~ 10% error level is comparable to the intrinsic 
uncertainty in the numerical models arising from our im- 
perfect knowledge of input parameters such as the H2 
formation rate coefficient and the dust cross section for 
LW photons. 

In the region % ^ 1, the two definitions begin to di- 
verge. Our analytic model still tracks the integral defini- 
tion reasonably well, reaching a maximum error of 40% at 
X = 0.078, the lowest x numerical run. In contrast, our 
model becomes increasingly inaccurate in estimating the 
location of the 50% molecular point as x ^ 0. This is not 
a surprising result. For % <C 1, we found in Paper I and 
we find again here that the transition between the atomic 
and molecular regions ceases to be sharp, and a major- 
ity of the atomic column is located in regions where the 
molecular fraction is > 50%. One can no longer identify 
a single meaningful location for the atomic-to-molecular 
transition. In this case our analytic model remains rea- 
sonably accurate in estimating the total atomic column, 
but, since we do not track the atomic fraction through 
the transition region, it does not provide a good esti- 
mate of where the atomic fraction reaches some specified 
value. This is analogous to the case of a Stromgrcn anal- 
ysis for ionizing radiation applied to a situation where 
the ionization front is not sharp (e.g. for ionization by x- 
rays rather the EUV photons) . In such a case there is no 
well-defined ionization front, and the Stromgren analysis 
does not give a good estimate for where the ionization 
fraction reaches 50% or some other specified value. How- 
ever, it does still provide a good estimate for the total 
emission measure, much as our method still gives a rea- 
sonably good estimate for the total H I column. 

Wc emphasize again that in Paper II we showed that 
X ~ 2 — 3 for any cloud, GMC or diffuse, in a region where 
two-phase equilibrium of the HI gas prevails. Values of 
X <C 1, where our method fails, can only be found in 
regions where the pressure is cither too low (e.g. well off 
the galactic plane, or in the far outskirts of galactic disks) 
or too high (e.g. in the nuclei of starburst galaxies) for 
two-phase equilibrium to be possible. 

4.2. Spherical Clouds 

Our algorithm for the spherical case with dust is a 
slightly more complex version of the ones described in 
§ 13.21 and § 14.11 We proceed as follows: 

1. Select values for r/j and Xfja- Our algorithm will 
determine the corresponding x- 

2. Pick a trial value of the parameter Xf. (The value 
of is fixed to 4/3, as discussed in § 13.11 ) The 
choice oi Xf determines the Eddington factor /Edd 
via equation (|5ip . 

3. Given /Edd, determine the function E* by integrat- 
ing the ordinary differential equations (|42)) and (|43p 
from X — to X = 1, using the boundary condi- 
tion that E*{xB.2) = 0. 

4. Check whether the solution obeys the integral con- 
straint equation (|54p to within some specified tol- 




FlG. 4. — Contours of xi, 

no 



as a function of and x- Solid 
lines show the numerical result, dotted lines show contours com- 
puted using the analytic approximation given by equation II82I I. 



The hatched region indicates whore x|j 



: 0. 



erance. If not, repeat steps 2-4 using Newton- 
Raphson iteration to search for a value of x/ 
that minimizes the difference between the left- and 
right-hand sides of equation (|54)) . 

5. If the solution does satisfy the constraint equation 
(|54p . determine x from E*{1) and r^? using equa- 
tion (I36l). 



We use this procedure to calculate values of x on a 
grid of values running from = 0.01 — 100 and X2 = 
0.01 — 0.99, and we use this table to construct contours 
of the molecular volume X2 in terms of x and r^j. The 
results are shown in Figurc21 As illustrated in the figure, 
the numerical result is well-approximated by 



t3 



1 



1 + 0.0712^, 



2.8 ' 



where 



ln(l-h0.6x + 0.01x^) 



0.6x 
71 fG'o 



ln(l -F0.6x + 0.01x^ 



0.6x 



(82) 



(83) 



(34) 



and where Nji,2i — uhR/ {10^^ cm ^) is the normalized 
column density in the cloud; recall that ^ = x/'TR- Note 
that the factor O.Olx^ is only important to fit the con- 
tours when X S> 1; and for x 1- the regime relevant 
to real molecular clouds, it may be omitted. Also note 
that equation (|82p is the same as equation (|70p for the 
dust-free spherical case, and the only difference between 
the two cases is that allows for a finite dust optical 
depth. In the limit ad (so that X ~* 0), we have 
£,d ^ S. and this result for xn^ approaches that found in 
i )3. 21 above for dust-free spheres. 

Altogether, we have results that cover most of the pa- 
rameter space in Figure [4] except for the region in and 
immediately adjacent to the fully atomic region. The 
results for a dusty sphere obtained in this section apply 
for X ~ 1 and th ^ 1 (the latter requirement is needed 
to ensure that the cloud is large enough for the dust to 
attenuate the radiation). For the region in which x < 1, 
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the results from WS.2\ are valid provided ^ — x/tr ^ 
but this is most of the region not adjacent to the fully 
atomic region. According to Paper II, the value of x in 
galaxies is typically ~ 2 — 3, depending on metallicity 
(eq. [HI). 

5. MOLECULAR MASS FRACTIONS IN 
ATOMIC-MOLECULAR COMPLEXES 

Thus far we have calculated molecular fractions under 
the assumption of uniform-density gas. However, it is 
convenient to generalize to the case where the atomic 
and molecular regions have different densities, since this 
is a more realistic representation of an atomic-molecular 
complex. Given a method to calculate function 
of and Tij, we solved the problem of determining the 
molecular mass fraction in such a cloud in Paper II. In 
that paper we used the approximate form for a;H2(Xi'''fl) 
determined in Paper I; in this section we apply the same 
approach to calculate the mass fraction for our new, im- 
proved determination of xu^ [x, tr). 

As discussed in Paper II, in a cloud where the atomic 
and molecular gas have differing densities, determining 
the molecular mass fraction requires an additional con- 
straint beyond the one imposed by radiative and chemical 
equilibrium. The most natural constraint is to require 
that the atomic and molecular gas be in pressure bal- 
ance. In an atomic-molecular complex of total (atomic 
plus molecular) surface density Scomp this implies that 



Tc = tr[1 



Omol 



l)x- 



where 



(85) 
(86) 



</'moi ~ 10 is the ratio of molecular to atomic gas vol- 
ume densities, and /xh is the mean mass per H nucleus. 
Given either an analytic approximation for X}i.^(xt'''r) 
(e.g. equation [5^ . or a numerical algorithm to calculate 
^H^iX: tr) (e.g. § SHI), equation ([85)l implicitly deter- 
mines the effective optical depth as a function of x- 
It is then simple to calculate and (x, ''i?) for a 
given value of x, and then to compute the corresponding 
molecular mass fraction 



/h2, 



^mol2;H2 



comp 



(0n 



1)^2 



(87) 



Alternately, one can obtain an approximate solution 
to equation (|85p analytically. Equations and (|87|) 
imply the simple relation 



rail - /h2, 



comp , 



= tr{1 



Solving this for /H2.comp and substituting in our approx- 
imate form for xn^ from equation (|82p gives 



/h 



2 , comp 



= 1 



tr 



1- 



where we have defined 



Tc V4/ 1 + 0.0712^^ 
3\ s 



2.8 



4 J l + 0.0712Cj ' 



TR ln(l+0.6x + 0.01x2) 
* = — ^d. = 



(89) 
(90) 

(91) 




100 

S„,„p (Mo pc-^) 

Fig. 5.- fn^ 

equation l|85| l. using the algorithm described in ? 14.21 to compute 
xii^ (x, ''"fi)- Dashed lines show the analytic approximation given by 
equation II93I I. Pairs of lines correspond to different metallicities, 
ranging from log Z' = —2 {furthest right) to log Z' = I {furthest 
left), in intervals of logZ' = 0.5. 

Note that s is defined solely in terms of given param- 
eters. The remaining unknown on the right-hand side 
of equation ([50]) is 0.0712^^-®, but since this is a small 
correction in the denominator, we can approximate it. 
Re-arranging equation (|85p . we find 



(I) 



(0„ 



(92) 



1+0. 0712^2. > 



and experimentation with numerical solutions to this 
equation shows that 0.0712^^'^ « 0.25s for (pmoi ~ 10. 
Substituting this into equation (j90p . we arrive at our ap- 
proximate solution: 



/h 



2 , comp 



1 + 0.25s 



(93) 



We apply this only for s < 2; for s > 2, we have 

f¥i2 , comp . 

Following Paper II, we can estimate x for a two- 
phase atomic medium using equation (jl3p . Similarly, 
we can write Tc in terms of the complex surface den- 
sity and the metallicity as Tc = 0.066I]comp,o^': where 
Scomp.o = Scomp/(l-^0 pc"^) and Z' is the metallic- 
ity relative to Solar. In writing this relationship, we 
have assumed the ad (x TZ (x Z' . With these esti- 
mate for X and Tc, we obtain the relationship between 
II2 mass fraction and Scomp shown in Figure [5] The 
difference between the approximate solution given by 
equation (1931) and the numerical solution is less than 



5% for all Scomp and metallicities between \ogZ' = 
—2 and logZ' = 1, where the difference between the 
approximate and numerical solutions is defined to be 



I/H2, 



approx 



./h2 



' max(/H2 1 - /H2,num)- 



6. CONCLUSIONS 



Using a Stromgren-type analysis, we have evaluated 
the thickness of the H I layer around molecular clouds for 
both planar and spherical geometry, and both without 
and with dust. Wc have assumed that the density of the 
HI layer is constant. In contrast with most treatments 
of PDRs, we have assumed that the radiation incident 
upon the cloud is isotropic. The results of our analysis 
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are based on several approximations. The major physi- 
cal approximation is that the gas in the layer is entirely 
atomic, /i = 1. We also assumed that the metallicity 
Z ^ Zq, so that the individual absorption lines in the 
Lyman- Werner bands overlap and the dissociati ng radi- 
ation covers most of the band 91.2-110.8 nm (iSternberel 
|1988() . The principal methodological approximation is 
the treatment of the variable Eddington factor, /Edd, as 
described in ^2.2\ In the absence of dust, f^dd for a 
spherical cloud has the unusual property that it is gen- 
erally less than the isotropic value of i. In addition, we 
assumed that the opacity associated with the flux was 
related to that associated with the number density by a 
constant, 0^. In the one case in which we co uld deter- 
mine as an eigenvalue of the solution fi i3.ip . we found 
that this was a good approximation. Finally, in calcu- 
lating the optical depth needed to obtain a constraint 
equation, we replaced the optical depth for individual 
rays, t/, with the angle-averaged value, te- The accu- 
racy of our approximations for the planar case has been 
tested by comparison with the results of the Meudon 
PDR code, which were kindly supplied by M. G. Wolfire. 
The analytic fit to our results (eq. [77]) typically agrees 
with the results of the PDR code to within 10%, with 
the maximum error of about 30%. 

Our results for the scaling behavior of the solutions 



are consistent with those of Paper 1. In the planar, dust- 
free case, the solution depends only on the characteristic 
Stromgren length, ^ch « Eq/u^TZ. When dust is in- 
cluded, the solution depends only on the dust optical 
depth associated with this length scale, x = '^HO'd^ch 
Gq/hh. In the spherical case, an additional parame- 
ter enters, the dust opacity associated with the cloud 
radius: th = niiadR- In the dust-free case, the so- 
lution depends only on the ratio of these two param- 
eters, ^ = x/tr = 4h/-R oc GyiriYiNR). We have 
presented simple analytic approximations for our results 
that should be useful for applications. 
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APPENDIX 
PRIMORDIAL FORMATION OF H2 
In the simplest case, H2 formation in low-density primordial gas is governed by just two reactions, the formation of 

H + c"^H"+lw, (Al) 

and the formation of II2 , 

H-+H^H2+e-. (A2) 

The reaction r ates for number of important reactions in primordial gas have been summarized recently by 
iGlover fc Abell (|2008[ ). The rate coefficient they cite for formation is within 15% of the simple expression 

fc_ = 1.83 X 10-18^0-88 cm3 s-i (A3) 

given by iHutchini (|1976f ) . They indicate that the rate coefficient they give for formation of H2 , 

fc2 = 1.3 X 10"^ cm3 s~\ (A4) 

is uncertain by at least a factor 2. If these two reactions are the only ones governing the abundance of H~ — in 
particular, if photodctachment by infrared radiation is negligible — then the equilibrium abundance of is 

n_ = ^1=-^ ne. (A5) 

The formation rate coefficient TZ- that is analogous to the rate coefficient TZ for II2 formation on dust grains, 

n^TZ- = k2?i-nY[, (A6) 

is then given by 

7^_ =8x lO-iV-s^a^'^ cm^s-i, (A7) 

where the normalized ionization, X p. -3 = n^/ (lO^^nn), is typical of regions of primordial star formation (e.g., 
lAbel et al.l[200l IBromm et al.l[200l .' 

While this analysis shows that in principle our formalism can apply to H2 formation in primordial clouds via the 
H" process, it is not clear whether it can do so in practice. For temperatures ^ 10^ K, the H~ rate coefficient is 
about 30 times less than the rate coefficient for H2 formation in the Milky Way. As a result, the time scale for II2 
formation is very long unless the density is high, but if the density is high it is likely that the ionization will be less 
than the assumed value of 10"'^ and that the formation-rate coefficient will be corr espondingly reduced. For example, 
in their calculations of radiative feedback by the first stars, lAhn fc Shapirol (|2007l ) found that the II2 concentrations 
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were generally less than about 1%, whereas our analysis determines the condition for the concentration to approach 
100%. 

At higher densities, H2 can form via three-body collisions at a rate 

niUsb = n^k^b, (A8) 

so that the rate coef hcient T^^ft is proport ional to the density. The three-body formation rate is uncertain by two 
orders of magnitude (jGlover fc Abell [2008( 1: the geometric mean of the high and low rates they cite is k^b = 2.4 x 
^Q_28y-i.27 pj^6 g-i fQj. > 3QQ corresponding to Usb = 3.7 x 10-^^{nii/10^° cm-^)T;^^ '^'^ cm^ s'^. Since 
three-body formation is important only at high densities, high UV fluxes, such as those in circumstellar disks or in 
close binaries, are required to create significant photodissociation zones. Such calculations have yet to be carried out. 
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